% Runs simluation of Big G model
% This version: September 2023
% Generates figures 5 and 6
% Calls mod file: modelG.mod

clear all; close all; clc

addpath c:\dynare\5.4\matlab

variables = char('federalpurchases','pi_c','i','y','g_output', 'c');
var_names = {'Federal purchases (percent)','Inflation (percent)','Interest Rate (p.p.)','Output (percent)','Spendingshare','Consumption (percent)'};



n_var = size(variables,1);


timemax = 21;
time    = 0:timemax-1;

dynare model_htm noclearall             % baseline model


options_.nomoments=1;
options_.nofunctions=1;
options_.irf=timemax;
options_.nograph=1;
options_.noprint=1;
IRF1=NaN(n_var,timemax);
IRF2=NaN(n_var,timemax);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% baseline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% new version August 2023
set_param_value('gy',.187);    % output share of total gov spending
set_param_value('rho_g1',.65); % persistence shock 1
set_param_value('rho_g2',.73); % persistence shock 2
set_param_value('std1',13.5);    % std shock 1 (percent)
set_param_value('std2',13.1);    % std shock 2 (percent)
set_param_value('alpha1',.78); % calvo sec 1
set_param_value('alpha2',.89); % calvo sec 2
set_param_value('gamma',.26);  % fraction of total spending going to sector 1 (interpolated form federal purchases)
set_param_value('n',0.65);     % size sector 1 (value added)

old = 0; % Version September 2022
if old 
    gamma = .56;
    gy = .2;
    set_param_value('gy',.2);  
    set_param_value('rho_g1',.85); 
    set_param_value('rho_g2',.85); 
    set_param_value('std1',1/gamma/gy/.16);  
    set_param_value('std2',1/(1-gamma)/gy/.16);  
    set_param_value('alpha1',.78);
    set_param_value('alpha2',.9);
    set_param_value('gamma',.56);
    set_param_value('n',0.72);  
end



%omega = (n - gamma*(1-zeta))/zeta;

[info, oo_] = stoch_simul(M_, options_, oo_, var_list_);
if info==0
    for varindex = 1:n_var
        IRF1_base(varindex,:) =  oo_.irfs.([ deblank(variables(varindex,:)) '_eps_g1']);
        IRF2_base(varindex,:) =  oo_.irfs.([ deblank(variables(varindex,:)) '_eps_g2']);
    end
else
    fprintf('Could not solve model for grid point %5.4f',ngrid(ijk))
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% symmetric model w/ average calvo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


averagecalvo = .5*.78+.5*.89;
set_param_value('gamma',.5);
set_param_value('n',.5);
set_param_value('alpha1',averagecalvo);
set_param_value('alpha2',averagecalvo);

if old 
    averagecalvo = n*.78+(1-n)*.90;
    set_param_value('gamma',.5);
    set_param_value('n',.5);
    set_param_value('alpha1',averagecalvo);
    set_param_value('alpha2',averagecalvo);
    gamma = 0.5;
    set_param_value('std1',1/gamma/gy/.16);   %.12
    set_param_value('std2',1/(1-gamma)/gy/.16);   %.24
end


[info, oo_] = stoch_simul(M_, options_, oo_, var_list_);
if info==0
    for varindex = 1:n_var
        IRF1_sym(varindex,:) = oo_.irfs.([ deblank(variables(varindex,:)) '_eps_g1']);
        IRF2_sym(varindex,:) =  oo_.irfs.([ deblank(variables(varindex,:)) '_eps_g2']);
    end
else
    fprintf('Could not solve model for grid point %5.4f',ngrid(ijk))
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure
for ijk = 1:5
    subplot(5,2,(ijk-1)*2 +1)
    plot(time,IRF1_base(ijk,:),'b',time,IRF1_sym(ijk,:),'r--','LineWidth',1.5); hold on;
    plot(time,zeros(1,timemax),'k--')
    title(deblank(var_names(ijk)) )
    ylabel(deblank(var_names(ijk)) )
  %   ylim([ymin(varindex) ymax(varindex)])
    xlim([0 timemax-1])
    if ijk == 1
        legend('Baseline','Symmetric model')
        title('Shock sector 1')
    end


    subplot(5,2,(ijk-1)*2+2)
    plot(time,IRF2_base(ijk,:),'b',time,IRF2_sym(ijk,:),'r--','LineWidth',1.5); hold on;
    plot(time,zeros(1,timemax),'k--')
   %  ylim([ymin(varindex) ymax(varindex)])
    xlim([0 timemax-1])
    if ijk == 1
        title('Shock sector 2')
    end

end



for ijk  =1:timemax
  discountvector(ijk) = 0.997^(ijk-1); 
end
discountvector = discountvector';

for ijk = 1:timemax
    pvmult1(ijk) = IRF1_base(find(strcmp(var_names,'Output (percent)')),1:ijk)*discountvector(1:ijk)./(IRF1_base(find(strcmp(var_names,'Spendingshare')),1:ijk)*discountvector(1:ijk));
    pvmult2(ijk) = IRF2_base(find(strcmp(var_names,'Output (percent)')),1:ijk)*discountvector(1:ijk)./(IRF2_base(find(strcmp(var_names,'Spendingshare')),1:ijk)*discountvector(1:ijk));
end


figure
plot(time,pvmult1,'b-',time,pvmult2,'r:','LineWidth',2);
legend('Sec 1 shock','Sec 2 shock','Interpreter','latex'); legend('BoxOff')
%ylim([0.35 .8])
%xlim([0 24])
grid on
xlabel('Months','Interpreter','latex')
set(gca,'TickLabelInterpreter','latex')
set(gca,'FontSize',30)
fig = gcf;
set(fig,'PaperOrientation','landscape');
print(gcf, ['../../figures/pvmult_htm.pdf'], '-bestfit', '-dpdf') 



figure
plot(time,IRF1_base(find(strcmp(var_names,'Spendingshare')),:),'k',time,IRF2_base(find(strcmp(var_names,'Spendingshare')),:),'k:',time,IRF1_base(find(strcmp(var_names,'Output (percent)')),:),'b-',time,IRF2_base(find(strcmp(var_names,'Output (percent)')),:),'r:','LineWidth',2); hold on
legend('Gov.spending Sec 1 shock','Gov.spending Sec 2 shock','Output Sec 1 shock','Output Sec 2 shock','Interpreter','latex'); legend('BoxOff')
ylabel('Percent of output','Interpreter','latex','FontWeight','normal');%,'FontSize',12,'FontWeight','normal');
%ylim([0 .35])
%xlim([0 24])
 grid on
xlabel('Months','Interpreter','latex')
set(gca,'TickLabelInterpreter','latex')
set(gca,'FontSize',30)
fig = gcf;
set(fig,'PaperOrientation','landscape');
print(gcf, ['../../figures/pvmult_components_htm.pdf'], '-bestfit', '-dpdf') 

%%
%% print each panel as a figure for paper
%%



ymax = [ 14 .1 .15 .35 .4 .1];
ymin = [ 0 -.0002  0 0  0 -.1];


for varindex = 1:n_var
    figure

        plot(time,IRF1_base(varindex,:),'b','LineWidth',5); hold on;
        plot(time,IRF1_sym(varindex,:),'r--','LineWidth',1); hold on;

    
    plot(time,zeros(1,timemax),'k:')
    ylim([ymin(varindex) ymax(varindex)])
    xlim([0 timemax-1])
    if varindex == 2
        legend('Baseline','Symmetric model','Interpreter','latex')
        legend BoxOff

     end
    grid on
    ylabel(var_names(varindex),'Interpreter','latex','FontSize',80,'FontWeight','normal');
    if varindex == 4
        xlabel('Months','Interpreter','latex')

    end
    
    set(gca,'TickLabelInterpreter','latex')
    set(gca,'FontSize',30)
    fig = gcf;
    set(fig,'PaperOrientation','landscape');
    print(gcf, ['../../figures/sec1_htm_' num2str(varindex) '.pdf'], '-bestfit', '-dpdf') 

    figure 

    plot(time,IRF2_base(varindex,:),'b','LineWidth',5); hold on;
    plot(time,IRF2_sym(varindex,:),'r--','LineWidth',1.5); hold on;
    
    plot(time,zeros(1,timemax),'k:')
    ylim([ymin(varindex) ymax(varindex)])
    xlim([0 timemax-1])
    grid on
     if varindex == 4
         xlabel('Months','Interpreter','latex'),
     
     end
    set(gca,'TickLabelInterpreter','latex')
    set(gca,'FontSize',30)
    fig = gcf;
    set(fig,'PaperOrientation','landscape');
    print(gcf, ['../../figures/sec2_htm_' num2str(varindex) '.pdf'], '-bestfit', '-dpdf') 
end
% close all